<?php
/**
 * Create by PhpStorm.
 * User: guokunsheng
 * Date: 2024/10/30
 * Time: 16:17
 * @author Seunk <85152470@qq.com>
 */

namespace app\common\lib\horoscope;

/**
 * 八字工具类
 * Class Util.
 */
class Util
{
    /**
     * 标准时间发出地经度(角度表示,东经为正西经为负),北京时间的经度为+120度0分
     */
    private $J = 120;

    /**
     * 默认纬度(角度表示,北纬为正南纬为负),这里是中国标准时间发出地(陕西省渭南市蒲城县)
     */
    private $W = 35;

    /**
     * 真太阳时模块,sn代表sin
     */
    private function sn($x)
    {
        return sin($x * 1.74532925199433E-02);
    }

    /**
     * 真太阳时模块,cn代表cosine
     */
    private function cn($x)
    {
        return cos($x * 1.74532925199433E-02);
    }

    /**
     * 真太阳时模块,返回小数部分(负数特殊) returns fractional part of a number
     */
    private function fpart($x)
    {
        $x = $x - floor($x);
        if ($x < 0)
        {
            $x = $x + 1;
        }
        return $x; //只取小数部份
    }

    /**
     * 真太阳时模块,只取整数部份
     */
    private function ipart($x)
    {
        if ($x == 0)
        {
            return 0;
        }
        return ($x / abs($x)) * floor(abs($x));
    }

    /**
     * 真太阳时模块,finds a parabola through three points and returns values of coordinates of extreme value (xe, ye) and zeros if any (z1, z2) assumes that the x values are -1, 0, +1
     */
    private function quad($ym, $y0, $yp)
    {
        $nz = 0;
        $A = 0.5 * ($ym + $yp) - $y0;
        $b = 0.5 * ($yp - $ym);
        $c = $y0;
        $xe = -$b / (2 * $A); //x coord of symmetry line
        $ye = ($A * $xe + $b) * $xe + $c; //extreme value for y in interval
        $dis = $b * $b - 4 * $A * $c; //discriminant
        $z1 = 0;
        $z2 = 0;
        if ($dis > 0)
        { //there are zeros
            $dx = 0.5 * sqrt($dis) / abs($A);
            $z1 = $xe - $dx;
            $z2 = $xe + $dx;
            if (abs($z1) <= 1)
            {
                $nz = $nz + 1;
            } //This zero is in interval
            if (abs($z2) <= 1)
            {
                $nz = $nz + 1;
            } //This zero is in interval
            if ($z1 < -1)
            {
                $z1 = $z2;
            }
        }
        return [$xe, $ye, $z1, $z2, $nz];
    }

    /**
     * 真太阳时模块,returns sine of the altitude of either the sun or the moon given the modified julian day of the UT
     * @param float $jd
     * @param float $J 经度,东经为正西经为负
     * @param float $W
     * @param int $LX 1月亮 2太阳日升日落 3太阳海上微光
     */
    private function sinalt($jd, $J, $W, $LX)
    {
        $instant = $jd - 2400001;

        $t = ($instant - 51544.5) / 36525; //减51544.5为相对2000年01月01日零点
        if ($LX == 1)
        {
            [$ra, $dec] = $this->moon($t);
        } else
        {
            [$ra, $dec] = $this->sun($t);
        }

        $mjd0 = $this->ipart($instant); //UT时间0点;returns the local sidereal time(计算观测地区的恒星时)开始
        $ut = ($instant - $mjd0) * 24;
        $t2 = ($mjd0 - 51544.5) / 36525;
        $gmst = 6.697374558 + 1.0027379093 * $ut;
        $gmst = $gmst + (8640184.812866 + (0.093104 - 0.0000062 * $t2) * $t2) * $t2 / 3600;
        $lmst = 24 * $this->fpart(($gmst + $J / 15) / 24); //结束

        $tau = 15 * ($lmst - $ra); //hour angle of object
        return $this->sn($W) * $this->sn($dec) + $this->cn($W) * $this->cn($dec) * $this->cn($tau);
    }

    /**
     * 真太阳时模块,关于太阳的,Returns RA and DEC of Sun to roughly 1 arcmin for few hundred years either side of J2000.0
     */
    private function sun($t)
    {
        $p2 = 2 * M_PI;
        $COSEPS = 0.91748;
        $SINEPS = 0.39778;
        $m = $p2 * $this->fpart(0.993133 + 99.997361 * $t); //Mean anomaly
        $dL = 6893 * sin($m) + 72 * sin(2 * $m); //Eq centre
        $L = $p2 * $this->fpart(0.7859453 + $m / $p2 + (6191.2 * $t + $dL) / 1296000);
        //convert to RA and DEC - ecliptic latitude of Sun taken as zero
        $sl = sin($L);
        $x = cos($L);
        $y = $COSEPS * $sl;
        $Z = $SINEPS * $sl;
        $rho = sqrt(1 - $Z * $Z);
        $dec = (360 / $p2) * atan($Z / $rho);
        $ra = (48 / $p2) * atan($y / ($x + $rho));
        if ($ra < 0)
        {
            $ra = $ra + 24;
        }
        return [$ra, $dec];
    }

    /**
     * 真太阳时模块,关于月球的,Returns RA and DEC of Moon to 5 arc min (ra) and 1 arc min (dec) for a few centuries either side of J2000.0
     * Predicts rise and set times to within minutes for about 500 years in past - TDT and UT time diference may become significant for long times
     */
    private function moon($t)
    {
        $p2 = 2 * M_PI;
        $ARC = 206264.8062;
        $COSEPS = 0.91748;
        $SINEPS = 0.39778;
        $L0 = $this->fpart(0.606433 + 1336.855225 * $t); //mean long Moon in revs
        $L = $p2 * $this->fpart(0.374897 + 1325.55241 * $t); //mean anomaly of Moon
        $LS = $p2 * $this->fpart(0.993133 + 99.997361 * $t); //mean anomaly of Sun
        $d = $p2 * $this->fpart(0.827361 + 1236.853086 * $t); //diff longitude sun and moon
        $F = $p2 * $this->fpart(0.259086 + 1342.227825 * $t); //mean arg latitude
        //longitude correction terms
        $dL = 22640 * sin($L) - 4586 * sin($L - 2 * $d);
        $dL = $dL + 2370 * sin(2 * $d) + 769 * sin(2 * $L);
        $dL = $dL - 668 * sin($LS) - 412 * sin(2 * $F);
        $dL = $dL - 212 * sin(2 * $L - 2 * $d) - 206 * sin($L + $LS - 2 * $d);
        $dL = $dL + 192 * sin($L + 2 * $d) - 165 * sin($LS - 2 * $d);
        $dL = $dL - 125 * sin($d) - 110 * sin($L + $LS);
        $dL = $dL + 148 * sin($L - $LS) - 55 * sin(2 * $F - 2 * $d);
        //latitude arguments
        $S = $F + ($dL + 412 * sin(2 * $F) + 541 * sin($LS)) / $ARC;
        $h = $F - 2 * $d;
        //latitude correction terms
        $N = -526 * sin($h) + 44 * sin($L + $h) - 31 * sin($h - $L) - 23 * sin($LS + $h);
        $N = $N + 11 * sin($h - $LS) - 25 * sin($F - 2 * $L) + 21 * sin($F - $L);
        $lmoon = $p2 * $this->fpart($L0 + $dL / 1296000); //Lat in rads
        $bmoon = (18520 * sin($S) + $N) / $ARC; //long in rads
        //convert to equatorial coords using a fixed ecliptic
        $CB = cos($bmoon);
        $x = $CB * cos($lmoon);
        $V = $CB * sin($lmoon);
        $C = sin($bmoon);
        $y = $COSEPS * $V - $SINEPS * $C;
        $Z = $SINEPS * $V + $COSEPS * $C;
        $rho = sqrt(1 - $Z * $Z);
        $dec = (360 / $p2) * atan($Z / $rho); //算出月球的视赤纬(apparent declination)
        $ra = (48 / $p2) * atan($y / ($x + $rho)); //算出月球的视赤经(apparent right ascension)
        if ($ra < 0)
        {
            $ra = $ra + 24;
        }
        return [$ra, $dec];
    }

    /**
     * 真太阳时模块,rise and set(升降计算) [升起时刻(真太阳时),落下时刻(真太阳时),真平太阳时差(仅类型2),升起时刻(标准时间,仅类型2),落下时刻(标准时间,仅类型2)]
     * @param float $jd
     * @param float $J 经度,东经为正西经为负
     * @param float $W
     * @param int $LX 类型:1月亮;2太阳日升日落;3太阳海上微光
     * @return array
     */
    private function risenset($jd, $J, $W, $LX)
    {
        $jd = floatval($jd);
        $J = +1 * floatval($J); //统一东经为正
        $W = +1 * floatval($W); //北纬为正,南纬为负

        $noon = round($jd) - $this->J / 360; //儒略日,中午12点,減去8小時時差

        $sinho = []; //太阳盘面几何中心与理想地平面之间的夹角
        $sinho[1] = $this->sn(8 / 60); //moonrise - average diameter used(月亮升降)
        $sinho[2] = $this->sn(-50 / 60); //sunrise - classic value for refraction(太阳升降)
        $sinho[3] = $this->sn(-12); //nautical twilight(海上微光)

        $rise = 0; //是否有升起动作
        $utrise = false; //升起的时间

        $sett = 0; //是否有落下动作
        $utset = false; //落下的时间

        $hour = 1;
        $zero2 = 0; //两小时内是否进行了升起和落下两个动作(极地附近有这种情况,如1999年12月25日,经度0,纬度67.43,当天的太阳只有8分钟-_-)

        $ym = $this->sinalt($noon + ($hour - 1) / 24, $J, $W, $LX) - $sinho[$LX]; //See STEP 1 and 2 of Web page description.
        if ($ym > 0)
        { //used later to classify non-risings 是否在地平线上方,用于判断极昼极夜
            $above = 1;
        } else
        {
            $above = 0;
        }

        do
        {
            //STEP 1 and STEP 3 of Web page description
            $y0 = $this->sinalt($noon + ($hour + 0) / 24, $J, $W, $LX) - $sinho[$LX];
            $yp = $this->sinalt($noon + ($hour + 1) / 24, $J, $W, $LX) - $sinho[$LX];
            //STEP 4 of web page description
            [$xe, $ye, $z1, $z2, $nz] = $this->quad($ym, $y0, $yp);
            switch ($nz)
            { //cases depend on values of discriminant - inner part of STEP 4
                case 0: //nothing  - go to next time slot
                    break;
                case 1: //simple rise / set event
                    if ($ym < 0)
                    { //must be a rising event
                        $utrise = $hour + $z1;
                        $rise = 1;
                    } else
                    { //must be setting
                        $utset = $hour + $z1;
                        $sett = 1;
                    }
                    break;
                case 2: //rises and sets within interval
                    if ($ye < 0)
                    { //minimum - so set then rise
                        $utrise = $hour + $z2;
                        $utset = $hour + $z1;
                    } else
                    { //maximum - so rise then set
                        $utrise = $hour + $z1;
                        $utset = $hour + $z2;
                    }
                    $rise = 1;
                    $sett = 1;
                    $zero2 = 1;
                    break;
            }
            $ym = $yp; //reuse the ordinate in the next interval
            $hour = $hour + 2;
        } while (!(($hour == 25) || ($rise * $sett == 1))); //STEP 5 of Web page description - have we finished for this object?

        if ($utset !== false)
        { //注意这里转成了真太阳时
            $utset = round($jd) - 0.5 + $utset / 24 - ($this->J - $J) * 4 / 60 / 24;
        }
        if ($utrise !== false)
        {
            $utrise = round($jd) - 0.5 + $utrise / 24 - ($this->J - $J) * 4 / 60 / 24;
        }

        $dt = 0; //地方平太阳时 减 真太阳时 的差值,即"真平太阳时差换算表",单位为天
        $tset = ($LX == 2) ? $utset : 0; //用于返回标准时间,关于月亮的必须先通过太阳升降获取到dt再转标准时间
        $trise = ($LX == 2) ? $utrise : 0;
        if (($LX == 2) && ($rise * $sett == 1))
        { //太阳相关,非极昼极夜且有升有落
            while ($tset < $trise)
            { //太阳先落下再升起,时区与经度不匹配的情况下会出现此种情况,加一天修正
                $tset += 1;
            }
            $dt = round($jd) - ($trise + ($tset - $trise) / 2); //单位为天.比较两者的中午12点(上午和下午是对称的)

            $tset = $tset - $dt + ($this->J - $J) * 4 / 60 / 24; //真太阳时转标准时间
            $trise = $trise - $dt + ($this->J - $J) * 4 / 60 / 24;
        }

        return [$utrise, $utset, $dt, $trise, $tset];
    }

    /**
     * 真太阳时模块,改编自 https://bieyu.com/ (月亮與太陽出没時間) 原理:用天文方法计算出太阳升起和落下时刻,中间则为当地正午(自创),与12点比较得到时差;与寿星万年历比较,两者相差在20秒内
     * @param float jd
     * @param float J 经度,东经为正西经为负,注意西经60度38分转换方式是: -60 + -1 * 38/60
     * @param float W 纬度,北纬为正南纬为负,太阳并不是严格从正东方升起,所以纬度也有影响,只是相对影响较小
     */
    private function zty($jd, $J, $W = null)
    {
        $jd = floatval($jd);
        $J = (is_null($J) === true) ? $this->J : floatval($J);
        $W = (is_null($J) === true) ? $this->W : floatval($W);

        [$utrise, $utset, $dt, $trise, $tset] = $this->risenset($jd, $J, $W, 2);

        return $jd - ($this->J - $J) * 4 / 60 / 24 + $dt; //转地方平太阳时+修正
    }

    /**
     * 將公历年月日時轉换爲儒略日历时间
     * @param int $yy (-1000-3000)
     * @param int $mm (1-12)
     * @param int $dd (1-31)
     * @param int $hh (0-23)
     * @param int $mt (0-59)
     * @param int $ss (0-59)
     * @return false|float
     */
    private function Jdays($yy, $mm, $dd, $hh = 12, $mt = 0, $ss = 0)
    {
        $yy = floatval($yy);
        $mm = floatval($mm);
        $dd = floatval($dd);
        $hh = floatval($hh);
        $mt = floatval($mt);
        $ss = floatval($ss);

        $yp = $yy + floor(($mm - 3) / 10);
        if (($yy > 1582) || ($yy == 1582 && $mm > 10) || ($yy == 1582 && $mm == 10 && $dd >= 15))
        {
            $init = 1721119.5;
            $jdy = floor($yp * 365.25) - floor($yp / 100) + floor($yp / 400);
        } else
        {
            if (($yy < 1582) || ($yy == 1582 && $mm < 10) || ($yy == 1582 && $mm == 10 && $dd <= 4))
            {
                $init = 1721117.5;
                $jdy = floor($yp * 365.25);
            } else
            { //不存在的时间
                return false;
            }
        }
        $mp = floor($mm + 9) % 12;
        $jdm = $mp * 30 + floor(($mp + 1) * 34 / 57);
        $jdd = $dd - 1;
        $hh = $hh + (($ss / 60) + $mt) / 60;
        $jdh = $hh / 24;
        return $jdy + $jdm + $jdd + $jdh + $init;
    }

    /**
     * 將儒略日轉换爲公历(即陽曆或格里曆) [年,月,日,时,分,秒]
     * @param float $jd
     * @return array
     */
    private function Jtime($jd)
    {
        $jd = floatval($jd);
        if ($jd >= 2299160.5)
        { //以1582年的10月15日0時(JD值2299160.5)為分界點,在這之前為儒略曆,之後為格里曆
            $y4h = 146097;
            $init = 1721119.5;
        } else
        {
            $y4h = 146100;
            $init = 1721117.5;
        }
        $jdr = floor($jd - $init);
        $yh = $y4h / 4;
        $cen = floor(($jdr + 0.75) / $yh);
        $d = floor($jdr + 0.75 - $cen * $yh);
        $ywl = 1461 / 4;
        $jy = floor(($d + 0.75) / $ywl);
        $d = floor($d + 0.75 - $ywl * $jy + 1);
        $ml = 153 / 5;
        $mp = floor(($d - 0.5) / $ml);
        $d = floor(($d - 0.5) - 30.6 * $mp + 1);
        $y = (100 * $cen) + $jy;
        $m = ($mp + 2) % 12 + 1;
        if ($m < 3)
        {
            $y = $y + 1;
        }
        $sd = floor(($jd + 0.5 - floor($jd + 0.5)) * 24 * 60 * 60 + 0.00005);
        $mt = floor($sd / 60);
        $ss = $sd % 60;
        $hh = floor($mt / 60);
        $mt = $mt % 60;
        $yy = floor($y);
        $mm = floor($m);
        $dd = floor($d);

        return [$yy, $mm, $dd, $hh, $mt, $ss];
    }

    /**
     *  计算地方真太阳时
     * @param $yy
     * @param $mm
     * @param $dd
     * @param $hh
     * @param $mt
     * @param $ss
     * @param float $J 经度
     * @param float $W 纬度
     * @return array
     */
    public function zhenTaiyangShi($yy, $mm, $dd, $hh, $mt, $ss, $J , $W)
    {
        $yy = (int)$yy;
        $mm = (int)$mm;
        $dd = (int)$dd;
        $hh = (int)$hh;
        $mt = (int)$mt;
        $ss = (int)$ss;

        $spcjd = $this->Jdays($yy, $mm, $dd, $hh, $mt, $ss); //special jd,这里依然是标准时间,即$this->J处的平太阳时

        //需要转地方真太阳时
//            $rt['pty'] = $spcjd - ($this->J - floatval($J)) * 4 / 60 / 24;
//            $rt['pty'] = $this->Jtime($rt['pty']); //地方平太阳时

        $spcjd = $this->zty($spcjd, $J, $W); //采用真太阳时排盘,这里有点疑问: 对应的廿四节气的计算是否也要转为真太阳时呢?
        $rt['zty'] = $this->Jtime($spcjd); //地方真太阳时

        return $this->Jtime($spcjd);
    }
}
